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Abstract 

This paper presents the numerical simulations of the Jet-A spray reacting flow in a single element 
lean direct injection (LDI) injector by using the National Combustion Code (NCC) with and without 
invoking the Eulerian scalar probability density function (PDF) method. The flow field is calculated by 
using the Reynolds averaged Navier-Stokes equations (RANS and URANS) with nonlinear turbulence 
models, and when the scalar PDF method is invoked, the energy and compositions or species mass 
fractions are calculated by solving the equation of an ensemble averaged density-weighted fine-grained 
probability density function that is referred to here as the averaged probability density function (APDF). 

A nonlinear model for closing the convection term of the scalar APDF equation is used in the presented 
simulations and will be briefly described. Detailed comparisons between the results and available 
experimental data are carried out. Some positive findings of invoking the Eulerian scalar PDF method in 
both improving the simulation quality and reducing the computing cost are observed. 

1.0 Introduction 

Many engineering applications of the computational fluid dynamics (CFD) for internal combustor 
flows need to accurately account for the turbulence-chemistry interaction to facilitate a higher fidelity 
analysis of the design. In the conventional CFD simulation methods that involve an averaging procedure 
(e.g., RANS and URANS) or filtering process (e.g., large eddy simulation LES), the turbulence-chemistry 
interaction manifests itself as an unclosed term that must be modeled. Various such models have been 
proposed. The simplest one is the so called “laminar chemistry” which simply ignores the effects of the 
fluctuations of turbulence on the chemical reactions. However, in the PDF simulation methods (no matter 
Fagrangian or Eulerian approach, joint or scalar PDF), the chemical reaction source term in the PDF 
equations is always in a closed form, therefore it can be directly calculated without any modeling (Ref. 1). 
This unique feature will stay when the PDF method is extended to its large eddy simulation (FES) 
approach (Ref. 2). 

In this paper, we present the preliminary numerical results of the complex spray reacting flow in a 
single element EDI injector. The simulation methods include both the conventional CFD methods 
(RANS, URANS) and an Eulerian scalar PDF method (Ref. 3), in which the velocity field is determined 
by the continuity and momentum equations of averaged Navier-Stokes equations. All simulations are 
done with the NCC code (Ref. 4) using the same numerical parameter setting and the same computational 
domain and grid resolution. In the RANS and URANS simulations, all the turbulent flux models, i.e., 
Reynolds stresses, heat and species fluxes, are nonlinear models (Ref. 5 and 6). The basic equations and 
models are described in Section 2.0. In the Eulerian scalar PDF method, we solve the transport equation 
of scalar APDF. It should be mentioned that the joint APDF equation for velocity and species is generally 


NAS A/TM— 20 1 2-2 1 7676 


1 



in a closed form (Ref. 7). However, for the marginal or scalar APDF equation, the convection term 
contains a conditional mean that needs to be modeled. Here a nonlinear model is introduced, which 
provides a diffusion process in the sample variable space. The brief description about the scalar APDF 
equation and model is given in Section 3.0. 

The liquid fuel is Jet-A, and C12H23 is adopted as its surrogate. At this stage of the simulations, we 
use a relatively simple spray injection model (i.e., prescribed droplet size and distribution) and a global 
five-species one-step kinetics for combustion chemistry (Ref. 8), so that we could concentrate more on 
the evaluation of different simulation approaches. 

The main results of steady simulations, including both with and without invoking an Eulerian scalar 
PDF method, are described in Section 4.1. They are compared side by side to examine how much effects 
are produced by the invoked scalar PDF method. The results for unsteady simulations are described in a 
similar manner in Section 4.2. 

One of the main objectives of this study is to search for a consistent and stable simulation tool for 
multiphase combustion flows. We require that this tool is not only to have the “physics” based models but 
is also able to produce the physically reasonable, numerically stable solution that can be sustained over a 
very long (infinite) time period. And also, this tool should be relatively economic. It is encouraging to 
observe from the present simulations that the adopted Eulerian scalar PDF method does show the 
potential for improving the simulation quality and also reducing the computing cost. 

2.0 Basic Equations for RANS and URANS Simulations 

2.1 Averaged Turbulent Variables 

In the case of compressible turbulent reacting flow, we often deal with two types of ensemble 
averaging: one with the density weighting, the other without the density weighting. The averaged 
turbulent variable without the density weighting is denoted by (|)(x,£) and is defined as 

< K *> 0 = lim T 7 X ' (*> ■ 0 ( ' 1 ) 

N-> oo TV 

n=l 

where (|) (w) is the n-th independent realization of turbulent variables, i.e., velocity components U\ n ^ , 

M 

density p (w) , pressure P (n \ m-th species mass fractions and internal energy e" n) = . N is the 

m = 1 

total number of independent realizations. The density-weighted averaged turbulent variable is denoted by 
§(x,t) and is defined as 


$(*»0 = -?r- ( 2 ) 

P 

These ensemble averaged variables (K*,0 , $(x,0 represent the turbulent mean variables that can be 
actually measured in experiments. They are deterministic in nature. They can vary with the time and 
space but usually contain relatively low frequencies and low wave numbers of the turbulent mean flow, 
when compared with the ones contained in the non-averaged turbulent variables (|>(x,f). 
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2.2 Basic Equation 

Applying the ensemble averaging, Equations (1) and (2), to the compressible Navier-Stokes equations, 
we will obtain: 


where 


dp , dp Uj Q 

dt dx , 


dpUi dp U i U j dp d 

1 — 1 

dt dxj dx t dxj 


2pv(S iJ --8 iJ S kk ) 


dpe G p U: e dq — — 

-¥—+ - ’ = — — + PS kk +2pv 
dt dx ; 0x, kk 


^ij^ij ^ S ii S kk 


\ 


+ Q 


j 


gpOm | dpU i O m _ a 

dt dx t dx t 


_p(m) d^m 
dX; , 




p = pr ±^=^ fptL, or P = 


m-l W m C v m = 1 


C v m = 1 




8X: 


i m = 1 


Sx 


(3) 

(4) 

(5) 

( 6 ) 

(V) 

( 8 ) 


In the above equations, k, v and T (w) are the molecular heat conductivity, kinematic viscosity and the m-th 
species diffusivity, they have the same dimension (i.e., velocity • length). It is commonly assumed that 
T (m) is the same for all species O m , and w m is the molecular weight. The h m , T are the enthalpy of species 
and the temperature, Q is the radiation rate, W m = p£ m is the chemical production rate of the m-th species, 
0 M + 1 represents the internal energy e , R is the universal gas constant. These equations are general; 
however, unlike the constant density flows, further approximations for the terms on their right hand side 
are required in order to complete the density weighted ensemble averaging process. One of such 
approximations leads to 


f 


2pv 


v 


S 5 -(v B 


A 




dUi dU, 


J 


dXj 


- + - 


dX: 


2 dud 

3 9 dx k 


(9) 


In which, we have basically neglected the variations of p and p during the averaging process, the value of 
p will be considered as a function of P, T, • • • . Other types of approximations are also possible, for 
example, 


2pv 


— SaSi 


ij^kk 


"dp Uj [ dpUj 
dx,- dx- 

V J 1 



dp ud 
d x k , 


( 10 ) 
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Similarly, 


_ dT ^ 3® 

q i =- pKC ^-2^ pr h > 


i m = 1 


8X: 


' ~ KC, 


dp T 

dX: 


M 

- X r(m) p 


m = 1 


do 

n i 
0X; 


(11) 


p(m) d0 m _ p(m) dp 0m 
dXj dXj 


( 12 ) 


Where v, k, c 0 , and r 1 '" 1 are considered as functions of P, T, • • • . Furthermore, we invoke the dissipation 
rate of turbulent kinetic energy, 


2pv 




- P 8 


Therefore, the ensemble averaged compressible Navier-Stokes equations can be written as 

dp [ dpUj _ Q 

dt dxj 


(13) 


(14) 


dpUj | dpUjUj _ dP | d 


dt 


dx • dx- dx • 

j 1 j 


r dpUi | dpUj 2 g 0p[/T 


r) V . 

v 


5x 




3x, 


oXi dx , 


3p® w 3p£/ ; ® w 


3^ 


- + 


3X 3X 


r(m) 3p®„ 


3x 


+ p5 w m = l,2,---,M 


i y 


p - = p,f v = £«f v, or 

|V C |V 

m=l w c v W =1 w 


- = p^|,0 m 0 


m^M+1 


C v m = 1 


(15) 


(16) 

(17) 

(18) 


<h = ~ Kc , 


0pT 

dX: 


M 


a ®» 

m = 1 


dX; 


(19) 


These equations are considered quite general, because i) they are exact if the flow becomes 
incompressible, ii) all the approximations made in Equations (15), (16) and (17) are related only to the 
molecular diffusion terms which are less important and even become negligibly small, i.e., they are of 
0(1 /Re) comparing with the convection terms on the left hand side for turbulent flows at high Reynolds 
numbers (see Ref 9 and 1). For the steady and unsteady simulations (i.e., RANS and URANS), Equation 
(14) to (17) are often used together with the further approximations for Equations (18) and (19): 


M 




or 


m = 1 


W w 


pRT I pRT 


M 


M 


1 _ ^p 0„, 1 _ ^p 0m 

M h w m ’ M £ W m 


( 20 ) 
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q.=-KC v ^— (21) 

OX i 

The momentum flux p £/,■£/• , the energy flux., and the species flux pU i O m are considered to be critically 

important and should be carefully modeled. Many models in the literature, from the simplest zero- 
equation model (Ref. 10) to the complex two-equation nonlinear models (Refs. 1 land 12), have been 

suggested in terms of turbulent stresses and turbulent scalar fluxes: x { - = p (£/,.£/. - Ufij) , 

0. = p (f/.O - f/.O) , where 0 represents the scalar quantities: e and O m . 

2.3 Nonlinear Closure Models for Turbulent Fluxes 

A general constitutive relationship between turbulent stresses x tj and the strain rate of mean 
turbulence^, Cl tj suggests (Ref. 13) 


k 3 , ~ ~ v 

P~^> -{Sik^kj - Qk S kj ) 

+2^5 P ~Sik^kj + ^ik^km^mj ~^kl 


S lm Q mk 8j3 + II s (S ij -d ij S kk /3)_, 


where, S y =(Uij +U U )/ 2, Cl tJ =(u iJ -Ujj)/ 2, II S =(V mm -S kl S lk )/l . The model 

coefficients A 3 and A 5 are constrained by the realizability condition and the rapid distortion theory 
limit. They are formulated as (see Ref. 12): 


Jl.O -4?c2 M \.6C^p— 

— A = J ^ ^ A — § 

''u h 9 3 , 2 5 5 7 4 _ ’ 

4.0 + 4 -t/* 0.5 + 1.5 ^qY 5 +QQ 

e e 2 s 3 4 


in which. 


^=V6coscp, cp = — arccos^V61F* j, IF 




\ tiA 0 ij°jk°ki 

»• w = a-f ■ 

(24) 

/ sa . 

(25) 


Note, here S- is a zero trace strain rate tensor. Similarly, the nonlinear model for scalar fluxes is 
formulated as (Ref. 14) 


0 . - 


<3p0 Q k> 
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Where denotes the turbulent diffusivity for the corresponding scalar quantity 0. It is often 
approximated by & r = v T / Pr e , and Pr e represents the turbulent Prandtl number or Schmidt number 
depending on whether the scalar quantity 0 is the energy e or the species . The turbulent eddy 
viscosity is defined as v r = • k 2 /s . The coefficients, C\ and c u are yet to be calibrated. In the current 

simulations they are set to be c x = c 2 = -0.24 

The turbulent kinetic energy and its dissipation rate k , 8 will be determined from the following model 
equations: 


d _ r d „ 7 0 

— p k H p Li- k = 

dt dx f 1 dx t 


(o + MrF k 

CX: 




P8 


d _ d _ „ d 

— P 8 H p it - 8 = 

dt dx ; dx ; 


dX: 


^ ~ e ^ pe 2 

- C o1 x lV v- C 0 


"ci ' 


"c2 " 


(27) 

(28) 


where C 8i and C e2 are the model coefficients. We have adopted the commonly used values of C 8 i = 1.45 
and C e2 = 1.92 in the present work while keeping in mind that they can be further constructed as functions 
of the local turbulence quantities (Ref. 15). 


3.0 Basic Equation for Scalar Averaged Probability Density Function 
(APDF) 

In this section, we will use the fine grained probability density function (FG-PDF) to define the 
averaged probability density function (APDF), then explore the relationship between APDF and the 
ensemble averaged turbulent scalar variables. This will provide the basis for establishing the transport 
equation for scalar APDF. 

3.1 APDF for Scalars 

3.1.1 Fine Grained Probability Density Function (FG-PDF) for Scalars (v|/; x,t) 

According to Pope’s definition (Ref. 1), the joint FG-PDF for turbulent velocity and scalars (i.e., 
compositions or species mass fractions, internal energy) can be written as 

3 M + 1 

/'(F,vi/;^0-5(^,0-F)8(O(x,0-v|/)-n 5 (V^0-^)II 8 ( (I) m( x ’0-VK m ) (29) 

1=1 m = 1 


Its marginal FG-PDF for scalars is 


M + 1 

/® O; X,t) = 8(®(*,0 - V) s ft 8(°m (*’0 “ Vm ) 

m= 1 


(30) 


where 8 denotes the delta function, U(x,t) is the turbulent (random) velocity vector (U u U 2 , C/ 3 ), O(x,0 is 
the turbulent (random) scalar array (0 1 , 0 2 , • • • , O m , O m+1 ) , for example, M species mass fractions and 
one internal energy O m+1 = e ; the xj denote the physical space variable (x u x 2 , x 3 ) and time t , V= (V u V 2 , 
V 3 ) and v| / = (\|/ 1? \|/ 2 ,-*-,\|/ M , \|/ M+1 ) are the sample space variables for U(x,t) and O (x,t), respectively. 
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3.1.2 Scalar APDF F 0 (v|/; x, t ) 

We define the following ensemble averaged density weighted fine-grained probability density 
function for scalars: 


x,t) = hm - v) 

n = 1 


( 31 ) 


Obviously, the scalar APDF, F 0 (v|/; x,t ) , is no longer a random quantity. And it satisfies the following 
“normalization” property: 

CO CO N , N 

J F (t ,(\\i-x,t)d\\i= | hm — ^p (w) (jc,^)8(o» ( " ) (jc,/) — \|/)^/\|/ = hm — yp ( " } (x,t) = p (32) 

-OO -00 v n=l v n = 1 


Note that here 



means that the integration is over the entire domain of the sample space \| /. 


3.1.3 Relationship Between Scalar APDF and Averaged Turbulent Variables 

With the definition of scalar APDF described in Equation (31), we can exactly deduce the averaged 
scalar turbulent variables that are defined in Equations (1) and (2). For example, 


+00 +°° f l N ^ 

f \|/F 0 (y; x,t)dy= f \\i hm — Yp ( yx,t)S(<I>(x,t)-x|/) <A|/ 

J - 00 J — GO N — >00 /V v ' 

V n=\ y 

V n=\ 




(33) 


1 N 

= lim — Y p ( ”l (x,t)<J>(x,t) 

n=\ 

= p® = p(x,t)®(x,t) 


Equation (33) denotes that the left hand side is an operation (o) that defines the density weighed 
turbulent mean variable p® : 


( 0 }=[ \|/F o (\|/;x,0<ty = p(*,0 <I> (x>0 

J —00 


(34) 


For a function W (®(x,t)) , it is easy to verify that 


(JE(®)) = J'VmF* (v|/; x,t) dy = p(x,t)W(<b(x,t)) 


(35) 


Furthermore, we may consider the derivatives V® as a new random quantity and legitimately write 


(V<D) = pVO 


(36) 
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However, because of the variable density, the “operation” ( ) does not have the differential commute 
property, i.e., 

(VO)^V(O) (37) 

because p VO =£ V(p<I>) . 

We can also write the density weighted mean (Uj 0 ( - \ for the joint velocity and scalar variables as 


(38) 


{Uj®i ) = J X v j Hh F u,<s> ( v > v; x,t)dVdyr = p(x,t)U j O i (x,t) 
where, F v 0 (K,v|/; x,t) is the joint APDF defined as 

F u<b (V,\ |/; x,t) = lim (x,t)5(c (n) ( x^-fWo^ (x,t ) - V ) 

N-> GO TV ; v ' 

n = 1 

which also satisfies the normalization property: 

00 00 1 ]\[ 

f f F U(J) (V,\\i;x,t)dVd\\i= lim — Y p (n) (x,i) = p 

J J ’ oo A/ “ 


(39) 


(40) 


n = 1 


3.1.4 Joint APDF, F c/ <1) (F,v|/; x,t) and Its Conditional APDF (F| v|/; x,t) 

From Equation (39), we may follow Reference 1 to define a “conditional” APDF on the condition 
O = \\r as 

F u<s> (V,W,x, t) 


F v ^{V | V ;x,0^-® 


^®(v; -M) 


(41) 


and the “conditional mean” as 


(U(x, 0|v} = |_7 v (F|v; x,t)dV = — j_ + ” F (F, y; x, f) JF 

00 X ^ t T ? J 00 


i i iv 

= ^Z p(n) (xM u(n) (*o ~ f K° ( ” ) (*o - v)f v 

= ^\im^^p {n \xd)U {n) §[& n \x,t)-y) 


(42) 


Then we have 


i ^ 

/7 = 1 


And the “complete” mean should be 
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J_ ((#(*, O|v)-O>)<ty 

= n_> «oi> w (at,<)-f) 5 (®« (*,,)- V )jvj v 

n = 1 

i N 

=^lim — ^p (n) (x,t)U {n) ( x,t ) = p(x,t)U(x,t ) 

n = 1 

Equation (43) can be extended to any other turbulent quantities, for example, VO, 5/(0): 
ton EgpW (V0) (n) 5(0 (fl) (** ) - x|/) = O • (vo|i|/) 

^Sp ( " ) 5 /" )8 (® (b) M - ■ v) = O • {s, (®)| V> = 0 ' ( v) 

«=§ 


(44) 


(45) 


where VO is viewed as a new random variable in addition to O. 


3.2 Transport Equation for Scalar APDF F 0 (y;x,t) 

We may derive the transport equation for scalar APDF from Equations (16) and (17) as follows: first, 
we write the terms on the left hand side of Equation (17) as 


op 0 /;? _ d 
dt dt 


= -\y m F <s> dy=\y m -^-dy 


(46) 


gp Up n 


SX: SX: 


{ f Vf y\> m 0,0 dVdy = 


—00 v l 

Then, the terms on the right hand side of Equation (17) can be written as 
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-00 -00 
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V dx i J 


dV d\]/ 


d\]f 


_d_ 

dX: 




r («) gpO 


3x 


* 7 


_d_ 

dX: 


T(m) jy]'Vm F * d 'V 


{ Vm 


—00 4 V 


p(m) 


8x 


(7v|/ 


i J 
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(47) 


(48) 


(49) 
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Where in Equations (48) and (49) we have applied the integration by parts and the following zero 
integration: 


00 

— nr> ' 


if a finite mean (Aj exists. 


( 50 ) 


Collecting all the integrand terms, and let the sum be zero (we refer this as the conservation condition, 
which is in general a sufficient but not necessary condition, Ref 16), we obtain 




dt 

or 


8X: 


dx, 


p(m) 


i V J 


8X: 




-[^•^(V)] , k = l,2,---,M + 1 


(51) 




dt 


dx , 


d d 

1 — 

dx Yk 3x ; - 


(m) llf 

dx- 




+ F 0 ^ -1,2,*-*,M + 1 


i J 


This equation also includes the equation of internal energy when S M+l (v|/) = 0 and other source terms in 
Equation (16) are neglected. Unlike the case of the joint APDF F v ^ , for the marginal F 0 , the 

convection term is not closed because of the conditional mean (U-| v|/^ . Then, this critically important 

term, corresponding to pLf-O in Equation (17), must be carefully modeled while the less important 
molecular diffusion term remains in the closed form. In addition, we noticed that the equally important 
chemistry source term p5 m (0) in Equation (17), which involves complex processes of turbulence- 
chemistry interaction, is closed in the scalar APDF equation, i.e., no need of modeling. This direct 
calculation of turbulence-chemistry interaction is one of the unique features of the PDF methodology. 


3.3 Nonlinear Model for Scalar Fluxes 


The convection term in Equation (51) contains the term • (U i | v|/) , which must be modeled. We 
may start from a nonlinear model, Equation (26), for the term pU i O m : 


pU i O m =pU i (S> m -r { 7 


(m) db;? 


dX: 


+ C 2 Qij ) 


dxj 


(52) 


This will lead to the following model by directly applying Equations (34) and (43): 




p(m) ^F^ 

T dx- 


+ - 


d Vk 


rp ) -(c l s 9 +c 2 n ij ) Wk 


8F® 
dx , 


or 


F <s> -(U i \y) = U i F (S) + 




p(m) 

T rk ~ 

dx t 


+ - 




<m) 


-(c\Sij +C 2 n,y)v|/ /f 


OF * 
dx- 


(53) 
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3.4 Summary 


With the model given by Equation (53), the scalar APDF equation for x,t) can be written as 


SF t | S(ViF») 

dt dx t 




+ 




8 





8F„ 


dx 


J J 


k = + 1 


(54) 


It can be verified that the scalar APDF Equation (54) can exactly deduce the ensemble averaged 
Equation (17). However, the model described by this equation is by no means unique. In addition, the 

variables ( Ui , Sy and ) are considered to be available during the solution procedure of the scalar 
APDF equation. Furthermore, in order to apply the available stochastic solution procedure built in the 
NCC code, we further simplify the model term in Equation (54) as follows: 


8 





dl% 

dx i j 


8 


1 p 
* 


\ 

J 


where 


\_ 

x 


F ( \, dx t 


r (m) 

i T 


~{ c l Sij +c 2 ~ V Sti Sij 


fix ■ 
UX J J 




(55) 


(56) 


Equation (56) is a crude approximation based on the dimensional argument for a time scale, which is 
responsible for the diffusion of scalar APDF in the space \| J k . We chose this time scale to be related to the 
rate of strain and rotation instead of its complex formulation. In order to prevent this time scale from 
being non-physically small during the simulation, we require that 

X > ^](v + v T )/e (57) 


because the right hand side of Equation (57) represents the smallest time scale of the simulated flow field. 


4.0 Numerical Simulations of Single Element LDI Injector 

The lean direct injection (LDI) injector is a liquid fuel injector developed to reduce aircraft emissions. 
Stable combustion is essentially completed within a short distance through rapid fuel and air mixing. This 
design also allows for many small fuel injectors integrated into modules facilitating different fuel staging 
strategies, such as the one shown in Figure 1. So far, experimental observations have not fully clarified 
the dynamics of the mixing and combustion processes occurring in these injectors, and numerical studies 
need to be conducted to achieve a better understanding of the underlying physics of the LDI combustor. 
Figure 2 shows the air s wirier and the convergent-divergent nozzle of the single element injector. Figure 3 
depicts the single element LDI combustor geometry and its computational domain. Five probes are 
dispatched along the centerline (Figure 4) to monitor the evolution of turbulent variables during the 
simulations. The numerical grid is formed using hexahedral elements and the total number of elements is 
about 862,000, which is a relatively coarse grid used in a previous RANS simulation (Ref. 17). 
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Figure 3. — Computational domain and grid spacing. 


Probe locations: 

1 . x = 0.00450m or -2.7 mm from dump face 

2. x = 0.01015m or 3 mm from dump face 

3. x = 0.01915m or 1 2 mm from dump face 

4. x = 0.13150m or 124 mm from dump face 

5. x = 0.1871 5m or 180 mm from dump face 



Figure 4. — Probes located on centerline. 


In this study, the liquid fuel is Jet-A, and C12H23 is adopted as its surrogate, the fuel is injected at the 
throat of the nozzle, mixing with the swirling air that comes from the air swirlers which consists of six 
helical, axial vanes with downstream vane angles of 60°. A prescribed droplet-size distribution spray model 
is used. In Section 4.1, we present the results of steady simulations. We first carry out a steady RANS 
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simulation with a “laminar chemistry” model, in which the generation rate of compositions is determined by 
the mean turbulent variables, which is commonly considered as a very crude assumption or model for the 
chemistry-turbulence interaction. This simulation has been carried out up to 900,000 iterations. Then we 
start a new simulation using the solution of the abovementioned simulation at the iteration of 460,000 as the 
initial condition and invoking an Eulerian scalar PDF method, in which the compositions and internal 
energy will be determined by the scalar APDF equation while the flow field is still determined by the 
continuity and momentum equations of the averaged Navier-Stokes equations. Since the production rate of 
composition is in a closed form in the scalar APDF equation, it can be calculated without modeling. This 
new simulation is carried out up to 760,000 iterations. In Section 4.2, we present the results of unsteady 
simulations URANS to consider the strong unsteadiness of EDI flow, using both with and without scalar 
PDF method. The stochastic numerical procedures for solving scalar APDF equation and fuel spray 
equations are described in Reference 3 and adopted in the present simulations. 


4.1 Results of RANS Simulations With and Without Invoking an Eulerian Scalar APDF 

In this section we present results of steady simulations with and without invoking the scalar APDF 
equation. Figure 5 shows the general pictures of simulated spray reacting flow in a single element FDI 
injector by both methods. Most of the results will be presented side by side for comparisons. The main 
results are presented in terms of 1) the variations of velocity components and temperature with respect to 
the iteration at five probes along the centerline (see Figure 4), which indicate how much the simulated 
reacting flow is developed, 2) the centerline distributions of temperature and axial velocity, 3) the 
temperature and velocity profiles along Z axis at several downstream locations, 4) the center recirculation 
zone visualized by isosurface of zero axial velocity and the contour plots of various turbulent quantities in 
the center x-z plane. These profiles, isosurface and contour plots will provide the information about flow 
and flame structures of the simulated spray reacting flow. Some available experimental data (Ref. 1 8) are 
also plotted for comparison with the calculated results from the simulations. 

4.1.1 Variation of Velocity and Temperature at Five Probes Along the Centerline 

Figure 6 and Figure 7 are the recorded iteration histories of velocity components and temperature at 
the five probes during the course of the steady simulations for the last 150,000 iterations. These figures 
indicate that the RANS simulation with the scalar APDF equation has reached its “steady” solution even 
before 600,000 iterations, the last 150,000 iterations are used for checking if the solution can be sustained 
or not. However, the conventional RANS without PDF method has the difficulty to reach its steady state 
and continues to vary around its “mean” value, which is clearly shown in the last 150,000 iterations. 
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LDI-spray-reacting flow 
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Figure 5. — Global pictures of spray reacting flow simulated with and without scalar PDF method. 
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Figure 6. — Comparison of variations of simulated velocity components at 5 probes. 
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Figure 6. — Concluded. 
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Figure 7. — Comparison of variations of simulated temperature at 5 probes. 
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Figure 7. — Concluded. 
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4.1.2 Centerline Distributions of Mean Temperature and Axial Velocity 

Figure 8 compares the centerline distribution of axial velocity with experimental data. The 
improvement in the nozzle region has been made by the PDF method. However, Figure 9 shows that the 
centerline distribution of temperature, comparing with the experimental data, is over predicted by the both 
methods in the region after the front face of dump combustor. 

4.1.3 Mean Temperature and Velocity Profiles Along Z Axis at Downstream Locations 

Temperature profiles along Z axis at ten downstream locations have been plotted. Here we only show 
four locations at x = 5, 10, 20 and 200 mm downstream of the front face of dump combustor. Comparing 
with the experimental data, Figure 10 clearly reveals the improvement made by the RAN simulation with 
the APDF equation. Similar situations are observed in the profiles of velocity components V and U along 
Z axis, shown in Figure 1 1 and Figure 12. 
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Figure 8. — Comparison of centerline distribution of axial velocity. 
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Figure 9. — Comparison of centerline distribution of temperature. 
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Figure 10. — Comparison of temperature distribution along Z axis at downstream locations. 
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Figure 10. — Concluded. 
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Figure 1 1 . — Comparison of V distribution along Z axis at downstream locations. 
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Figure 12. — Comparison of U distribution along Z axis at downstream locations. 
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Figure 13. — Comparison of center recirculation zone. 
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4.1.4 Mean Flow Structure and Contour Plots of Variables in the Center x-z Plane 

To reveal the simulated flow and flame structures, we have plotted the center recirculation zone 
visualized by the isosurface of zero axial velocity (see Figure 13), the contour plots of velocity 
components, temperature, turbulent kinetic energy, species mass fraction C12H23, 02 and C02 are 
shown in Figure 14, Figure 15, and Figure 16. From these figures we observed that all the structures 
simulated by RANS with APDF equation are more symmetric than the ones from the RANS without PDF 
method. This trend of symmetry seems more lined up with the experimental observation. 
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Figure 14. — Comparison of U,V,W contours in center plane. 
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Figure 15. — Comparison of temperature and k contours in center plane. 
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4.2 Results of URANS Simulations With and Without Invoking an Eulerian Scalar 
APDF 

In this section, we present the results of unsteady simulations with and without scalar APDF equation. 
Most results will be presented in a similar way to what we have done for the steady simulations in 
Section 4.1. URANS without the scalar APDF equation is carried out up to 83,000 time steps, the size of 
time step dt = 2.x 10 -6 , which is also used in the solvers for spray and scalar APDF equation for 
synchronizing the simulation. After the simulated flow is fully developed, the mean turbulent quantities 
are calculated by averaging over the last 20,000 time steps, which is about 2.44 “flow through times” 
(FTT) which is defined here as the time period that is needed for the flow to go from the inlet to the outlet 
of the computational domain based on the speed at the inlet. Then the simulation is switched to URANS 
with the scalar APDF equation and carried out up to 153,000 time steps. Again, the mean values of 
turbulent quantities are obtained by averaging over the last 20,000 time steps. Figure 17 depicts the global 
mean structures from abovementioned two unsteady simulations of spray reacting flow. 
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The main results to be presented are: 1) the convergence history of unsteady simulations, 2) the time 
histories of velocity components and temperature at four probes along the centerline, 3) the centerline 
distribution of temperature and axial velocity, 4) the profiles of temperature and velocity components 
along Z axis at several downstream locations compared with available experimental data, 5) the center 
recirculation zone and the contour plots of various mean turbulent quantities in the center x-z plane. 

4.2.1 Convergence History 

Figure 18 shows the time history of the number of subiteration incurred in unsteady simulations with 
and without invoking the scalar APDF equation over the last 20,000 time steps. It is interesting to see that 
the simulation with the scalar APDF equation converges much faster than the simulation without the 
scalar APDF equation. For each time step, the former converges after 21 to 33 subiterations, but the later 
often needs much more subiterations to converge. Comparing the wall time of computing, the simulation 
with APDF is about 31 hr with 256 processors for the 20,000 time steps; however, the simulation without 
APDF needs about 42 hr to finish 20,000 time steps of calculations. 
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Figure 17. — Global pictures of spray reacting flow simulated with and without scalar PDF method. 
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Figure 18. — Convergence histories of unsteady simulations. 
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4.2.2 Time Variation of Instantaneous Variables at Centerline Probes 


Time history of velocity components and temperature at four locations are recorded during the 
simulation. Here we list the time history at the probes 1, 2, 3, 5 (see Figure 19 and Figure 20), it can be 
seen that at the first three locations the flow appears to be strongly unsteady and fully developed. 
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Figure 19. — Time history of velocity components u, v and w at 4 probes. 
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Figure 19. — Concluded. 
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Figure 20. — Time history of temperature at 4 probes. 
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4.2.3 Centerline Distributions of Mean Temperature and Axial Velocity 

Centerline distributions of mean axial velocity and temperature are shown in Figure 20 and Figure 21. 
The effect of scalar APDF equation on the unsteady simulation is quite similar to the case of steady 
simulation described in Section 4.1.2. 

4.2.4 Mean Temperature and Velocity Profiles Along Z Axis at Downstream Locations 

Mean temperature profiles at downstream locations x = 5, 10, 20 and 200mm are compared and 
shown in Figure 22. These comparisons have revealed a positive effect of the APDF equation on the 
prediction of temperature. However, the effect of APDF equation on the velocity field, shown in Figure 
23 and Figure 24, is not as obvious in present unsteady simulations. 
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Figure 22. — Centerline distribution of mean temperature T. 
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Figure 23. — Comparison of temperature distribution along Z axis at x = 5, 10, 20, 200 mm. 
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Figure 23. — Concluded. 
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Figure 24. — Comparison of velocity component V along Z axis at x = 5, 15, 29, 38, 92mm. 
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Figure 24. — Concluded. 


4.2.5 Flow Structure and Contour Plots of Mean Variables in the Center x-z Plane 

To reveal the simulated mean flow and flame structures, we have plotted the center recirculation zone 
by the isosurface of zero axial velocity (see Figure 25), the contour plots of velocity components, 
temperature, turbulent kinetic energy and species mass fractions C12H23, 02 and C02 are shown in 
Figure 26, Figure 27, and Figure 28. From these figures we observe that the structures of all scalar fields 
simulated by URANS with the APDF equation are quite different from the ones by URANS without the 
APDF equation. Also, we notice that the level of turbulent kinetic energy is significantly reduced in the 
unsteady simulation with the APDF equation (a phenomenon often found in the simulations with the 
conventional or standard PDF method); however, it does not severely affect the global structure of 
turbulent reacting flow in our simulations due to, we believe, the use of the APDF equation, in which the 
time scale or frequency for the APDF diffusion in the sample space is not directly related to the turbulent 
kinetic energy, instead, it is determined by the strain and rotation rate of the turbulent mean flow (see 
Eq. (56)). 
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Figure 25. — Comparison of velocity component U along Z axis at x = 5, 15, 29, 38, 92mm. 
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Figure 25. — Concluded. 
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Figure 26. — Comparison of center recirculation zone. 
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Figure 27. — Comparison of mean velocities U, V, W contours in center plane. 
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Figure 28. — Comparison of mean temperature and k contours in center plane. 
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Figure 29. — Comparison of C12H23, 02 and C02 contours in center plane. 


5.0 Conclusions 

The numerical simulations of Jet-A spray reacting flow in a single element LDI injector have been 
carried out both with and without invoking an Eulerian scalar APDF method. The NCC code is used 
under the same numerical setting, same computational geometry and grid resolution. In addition, the same 
spray model and the same chemistry kinetics are used for all the steady and unsteady simulations. In this 
way, we hope to isolate and study only the effect of scalar PDF method (with the scalar APDF equation) 
on simulations, paying particular attention to the simulation quality and numerical performance (e.g., 
computational stability and efficiency). 

From both steady and unsteady simulations we have observed that there are noticeable improvements 
achieved on the temperature prediction in the region of strong turbulence by using the scalar PDF method. 
The distribution or structure of other scalar quantities appears to have significant differences between the 
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simulations using with and without the scalar APDF equation, and the scalar PDF method seems to 
produce more reasonable results. However, we need more reliable experimental data for a better 
evaluation. For the flow field, we also notice a big difference and improvement in the steady simulation 
with the PDF method, but not much in the unsteady simulation. 

As to the numerical performance, the observation is that both steady and unsteady simulations 
invoking an Eulerian scalar PDF method appear to be more stable and faster than the conventional RANS 
and URAMS simulations. 

Further evaluations of the scalar APDF equation and studies of stochastic numerical procedure for 
Eulerian scalar PDF method will be continued and extended to the large eddy simulation (LES) approach. 
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